set.seed(10551)

rm(list = ls())

library(nnet)
library(modelsummary)
library(tidyverse)

## Load MP data (including abstensions/no shows)

dat <- readRDS('data/mp_data.rds')  %>%
  filter(elecper != 8) %>%
  mutate(outcome_cat = case_when(vote_yes == 1 ~ 'yes',
                                 vote_abstain_explicit == 1 ~ 'abstain',
                                 vote_abstain_no_show == 1 ~ 'absent',
                                 vote_no == 1 ~ 'no'))

## Relevel categorical outcome

dat$outcome_cat <- relevel(as.factor(dat$outcome_cat), ref = 'no')
dat$outcome_cat <- as.factor(dat$outcome_cat)

## Define covariates 

experience_covars <- c('nsdap_member', 
                       'syn_in_gem_bin', 
                       'rel_cath',
                       'veteran_ww1',
                       'veteran_ww2', 
                       'capture_ww2',
                       'soviet_capture',
                       'exile_repressed_combined', 
                       'resistance_member')

covars <- c('year_birth', 'gender', 'mandate', 'dualcand', 'closeness_district_categ')

## Estimate multinomial models

m1 <-  paste('outcome_cat', '~', paste(experience_covars[2], collapse = "+"),' + elecper')
m2 <-  paste('outcome_cat', '~', paste(experience_covars[2], collapse = "+"),' + elecper + party_fe')
m3 <-  paste('outcome_cat', '~', paste(experience_covars[2], collapse = "+"),' + elecper + party_fe + state_id')
m4 <-  paste('outcome_cat', '~', paste(experience_covars, collapse = "+"),'+', paste(covars, collapse = "+"),' + elecper')
m5 <-  paste('outcome_cat', '~', paste(experience_covars, collapse = "+"),'+', paste(covars, collapse = "+"),' + elecper + party_fe')
m6 <-  paste('outcome_cat', '~', paste(experience_covars, collapse = "+"),'+', paste(covars, collapse = "+"),' + elecper + party_fe + state_id')

mlist <- list(m1, m2, m3, m4, m5, m6)

## Estimate 

m_est1 <- multinom(m1, data = dat, maxit = 1000) %>% modelsummary::get_estimates()
m_est2 <- multinom(m2, data = dat, maxit = 1000) %>% modelsummary::get_estimates()
m_est3 <- multinom(m3, data = dat, maxit = 1000) %>% modelsummary::get_estimates()
m_est4 <- multinom(m4, data = dat, maxit = 1000) %>% modelsummary::get_estimates()
m_est5 <- multinom(m5, data = dat, maxit = 1000) %>% modelsummary::get_estimates()
m_est6 <- multinom(m6, data = dat, maxit = 1000) %>% modelsummary::get_estimates()

## To latex Table

m_res <- bind_rows(m_est1, m_est2, m_est3, m_est4, m_est5, m_est6) %>%
  filter(term == 'syn_in_gem_bin' & response == 'yes') %>%
  dplyr::select(estimate, std.error, p.value) %>%
  mutate_all(~ round(., 4)) %>% 
  mutate(covariates = c(rep('No', 3), rep('Yes', 3)),
         fe = rep(c('Election period', 'Election period + Party', 'Election period + Party + State'), 2)) %>%
  kableExtra::kable(., 
                    col.names = c('Estimate', 'Std. error', 'P-value', 'Covariates', 'Fixed effects'),
                    format = "latex", booktabs = TRUE, caption = "Results from multinomial regression including absentees")


m_res

